data = readtable('../../data/main_VAR/UK_source_pre1946.xlsx', 'Sheet', 'Full_Data_for_VAR');
date = data.Year;
date0 = date;

startdatenum = 1729;
enddatenum = 1946;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data(find(date0 >= startdatenum - 1 & date0 <= enddatenum), :);

taxrevgdp = data.revtogdp(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.revtogdp(1); % government tax revenue to GDP ratio
spendgdp = data.spendingtogdp(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.spendingtogdp(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.realGdpGrowth(2:end); % real gdp growth (in logs)
inflation = data.inflation(2:end); % growth in GDP price deflator (in logs)
ynom1 = data.short_r(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.rate_10year(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

dy = data.dy(2:end);
pdm = log(1 ./ dy); % log price/dividend ratio

%% convenience yield

cy_short1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck5:ck45'); % 1873-1913
cy_short2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck57:ck63'); % 1925-1931

cy_short = [mean(cy_short1) * ones(1873 - startdatenum + 1, 1); cy_short1; mean(cy_short2) * ones(1925 - 1914, 1); ...
                cy_short2; mean(cy_short2) * ones(enddatenum - 1931, 1)] / 100;

cy_long1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck5:ck45'); % 1873-1913
cy_long2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck57:ck63'); % 1925-1931

cy_long = [mean(cy_long1) * ones(1873 - startdatenum + 1, 1); cy_long1; mean(cy_long2) * ones(1925 - 1914, 1); ...
               cy_long2; mean(cy_long2) * ones(enddatenum - 1931, 1)] / 100;

cy = cy_long;

ynom1 = ynom1 + cy_short(2:end);
ynom10 = ynom10 + cy_long(2:end);

ynom1 = log(ynom1 +1);
ynom10 = log(ynom10 +1);

divgrm = diff(log(data.dy) + data.stockindex);

gdebt = data.debttogdplevel;
% calibrate cy/gdp ratio
cy_data = mean([(cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2);
                (cy_long2) / 100 .* gdebt(find(date == 1925) + 1:find(date == 1931) + 1)]);
cy_data1 = mean((cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2));
cy_data2 = mean((cy_long2) / 100 .* gdebt(find(date == 1925) + 1:find(date == 1931) + 1));

seig_rev = [cy_data1 * ones(1873 - startdatenum + 1, 1); gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2) .* (cy_long1(1:end) / 100); ...
                ones(1925 - 1914, 1) * cy_data2; gdebt(find(date == 1925) + 1:find(date == 1931) + 1) .* cy_long2 / 100; cy_data2 * ones(enddatenum - 1931, 1)];

% adjust for convenience yield
taxrevgdp = taxrevgdp + seig_rev(2:end);
taxrevgdp0 = taxrevgdp0 + seig_rev(1, 1); % use gdebt in year 1729 as the value in 1728;

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1; % spread between 10-yr and 1-yr yields
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tts(1) = taxrevgdp(1);
gts(1) = spendgdp(1);

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

%% 
run ../../tools/run_var_estimate;

save(['MAT/res_annual_VAR_nodebt_cy.mat'], '-regexp', '^(?!(option|Globaloption)$).');
